A fluctuating lattice-Boltzmann method with improved Galilean invariance 
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In this paper we show that standard implementations of fluctuating Lattice Boltzmann methods 
do not obey Galilean invariance at a fundamental level. In trying to remedy this we are led to a 
novel kind of multi-relaxation time lattice Boltzmann methods where the collision matrix depends 
on the local velocity. This new method is conceptually elegant but numerically inefficient. With 
a small numerical trick, however, this method recovers nearly the original efficiency and allows 
the practical implementation of fluctuating lattice Boltzmann methods with significantly improved 
Galilean invariance. This will be important for applications of fluctuating lattice Boltzmann for 
non-equilibrium systems involving strong flow fields. 



I. INTRODUCTION 



Including noise in lattice Boltzmann simulations has 
been an active field of research in the last few years. 
It was pioneered by Ladd[l] who suggested to introduce 
noise on the non-conserved hydrodynamic modes, i.e. the 
stress modes. This approach works reasonably well in 
the hydrodynamic limit but for short length scales the 
fluctuations are underrepresented due to interaction with 
the non-hydrodynamic degrees of freedom which are typ- 
ically called the 'ghost'-modes. Adhikari et al. [2| recog- 
nized the necessity to include noise on all non-conserved 
degrees of freedom, including the non-physical 'ghost'- 
modes and Diinweg et al. [3] reformulated this approach 
to follow a detailed-balance condition description. All of 
these publications describe a fluctuating isothermal ideal 
gas. Just recently there was significant progress in ex- 
tending this concept to non- ideal equations of state (343 • 

The Adhikari implementation employs a multi- 
relaxation time (MRT) method similar to the one origi- 
nally introduced by d'Humieres Q except that the modes 
are orthogonal with respect to the Hermite norm. This 
allows for independent relaxation to the physically rele- 
vant moments. In particular it simplifies the construc- 
tion of a noise term that does not violate conservation 
laws while allowing for non-correlated noise on all other 
degrees of freedom. The derivation of the fluctuation- 
dissipation theorem in both, Adhikari's and Diinweg's 
approaches requires the MRT transforms to be orthog- 
onal with respect to a certain norm. In the case of a 
fluctuating ideal gas this norm depends on the equilib- 
rium distribution. However, in all previous publications 
the equilibrium distribution in this norm is taken only to 
zeroth order, i.e. only the weight factors in the equilib- 
rium distribution are used. The result is that the MRT 
orthogonality condition employed is identical to what is 
typically known as the Hermite norm Q. This approxi- 
mation, as we first discussed in Q and show later in this 
paper, formally introduces non-Galilean invariant terms. 
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We investigate here the effects of using this zeroth order 
approximation with respect to fluctuations in the context 
of non-zero flow speeds. The observed Galilean invari- 
ance violations suggest that this approximation may be 
inappropriate in some cases. To avoid this approximation 
we developed a novel kind of lattice Boltzmann method 
which includes the full second order expression which we 
expected to significantly reduce the Galilean invariance 
violations observed. Such a method necessarily has a lo- 
cal collision matrix that depends on the velocity at the 
respective lattice site. 

The paper is structured as follows: In section two we 
present a more detailed derivation based on Adhikari's 
noise implementation to show where the non-Galilean in- 
variant terms originate. We elaborate on the source of 
the orthogonality condition and the consequences of the 
zeroth order approximation and illustrate the impact on 
the MRT transforms. In section three we test the current 
literature standard for the example of a D2Q9 simula- 
tion. We measure the validity of two core assumptions 
of the derivation in the context of large flow speeds and 
find that Galilean invariance is indeed violated. Section 
four then discusses approaches to remedy the Galilean 
invariance violations. In particular we move away from 
the zeroth order orthogonality condition and attempt to 
introduce first and second order velocity terms of the 
equilibrium distribution. As a consequence we derive a 
lattice Boltzmann method for which the MRT transforms 
become locally velocity dependent. However, a simplistic 
implementation of this method is numerically inefficient. 
This inefficiency can be overcome by introducing look-up 
tables. The resulting LB scheme's computational cost is 
only slightly larger than that of the Hermite norm imple- 
mentation and Galilean invariance violations are signifi- 
cantly reduced. 



II. LATTICE BOLTZMANN SIMULATION OF A 
FLUCTUATING IDEAL GAS 

In order illustrate the origin of Galilean invariance vio- 
lations in fluctuating lattice Boltzmann implementations 
we present a short derivation of the fluctuating ideal gas 
in the Lattice Boltzmann context. The derivation pre- 
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sented is based on Adhikari et al.'s work [2] who first 
recognized the necessity to include noise on all non- 
conserved degrees of freedom. The derivation given in 
Adhikari et aVs original paper is not very detailed and 
we clarify some of the omitted steps of their derivation in 
this section. We put emphasis on a clear notation that 
separates the velocity space distibution functions fi and 
the moment space moments we call M a . 

The fluctuating lattice-Boltzmann equation is given by 

/i(x+ Vi,t + 1) = (1) 
/,•(*,<) + £A y [fj(x,t) - /?(x,i)] + fi(x,t), 

3 

where the fi are densities associated with the velocities 
Vi. The local equilibrium distribution depends on posi- 
tion and time through the local density p=Y,ifi an d ve- 
locity u = Y,i fiVi/p. The structure of the collision matrix 
Ay is discussed later in this section. This is the standard 
BGK lattice-Boltzmann equation with an added noise 
term £j(x, t). These noise terms must be chosen such 
that conserved quantities p, j, where j = £j/jV,j, a re 
not changed and a proper fluctuation dissipation theo- 
rem (FDT) is obeyed. How we obtain the latter while 
ensuring the former is outlined below. 

Throu gho ut this paper we use Qian's second order ex- 
pansion [1C| of the continuous Maxwell-Boltzmann dis- 
tribution as expression for the equilibrium distribution 



/°(p,M) 
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(2) 



This form is typically used for simulations of isothermal 
hydrodynamics. The extention to thermal hydrodynam- 
ics is conceptually straight forward. All references below 
to zeroth, first or second order terms in velocity of the 
equilibrium distribution are to be understood in terms 
powers of u of this expression. 

In order to gain independent access of conserved and 
non-conserved moments it is useful to shift from Boltz- 
mann type particle distributions /, to what is called 
generalized lattice-Boltzmann, moment space represen- 
tation, or multi relaxation time representation (MRT)[7, 
fl~0 | . One thus gains access to the hydrodynamically rele- 
vant moments directly. For this purpose a set of a forward 
transform from velocity space and its density functions 
fi to moment space and its so-called moments M a 



M°(x,i) = ^m J a / i (x,t). 

i 

and the corresponding back transform 
/,(x,t) = £<M»(x,t). 



(3) 



(4) 



must be chosen. While the original matrix elements m° 
and nf in were identical this is not necessary. But 
they need to follow the orthogonality conditions 



The particular choice of these transforms aims to gen- 
erate a simple form for the fluctuation dissipation theo- 
rem and is of key importance to the validity of the noise 
derivation and Galilean invariance or lack thereof. As 
such they differ from those in the publications introduc- 
ing the MKT formalism 0, HH • At least in the case of 
the ideal gas implementation it is convenient to choose 
the moments M a such that the representation of the col- 
lision matrix A in moment space is diagonal A afc = -^S ab . 
For practical purposes it is then useful to perform the 
collision in moment space. The fluctuating LBE Eq. ([T]) 
is then written as 



/i(x + u i5 i+l) - /i(x,t) = 

£ < h A° b [M b (x, f) - M fe '°(x, 0] + r 



(6) 



where £ a is the noise amplitude associated with moment 
M a and A is a random number chosen from a Gaus- 
sian distribution with a variance of one. The primary 
advantage here is that we gain independent access to the 
hydrodynamically relevant physical moments and we can 
choose the noise amplitudes £ a such that conservation 
laws are not violated, i.e. ^».«m«r«ed = 

Now we separate the fi in Eq. (JTJ) into their global 
mean values and a local fluctuating term 



fi = {fi)+Sfi 



(7) 



and we obtain 



(f i ) + Sf i (x + v i) t+l) = (f i ) + Sf i (y:,t) (8) 
+ lA ii [(/ j > + i/ j (x,t)-{/»}-^(x ) t)] 

3 

+ fc(x,t). 

Subtracting the (/,•) and assuming 

(fi} = f?(p ,u ), (9) 

where po and Uo are the equilibrium values of the density 
and the velocity, yields a LBE for the fluctuation part of 
the distribution 

<y/ i (x+« il t+i)= (io) 

Sf l (x,t) + Y l A lJ [*/ i (x J t)-*/9(x,i)]+&(x J t). 

3 

We can now Fourier transform in space and apply the 
moment space transform £i m° to obtain the moment 
space evolution equation in fc-space 

SM a (k,t + 1) = Y J Y J < e ' lkVtn \{ 5Mb ( k ^) + ( n ) 

i b 

ZEE AbCm 3 n 3 [5M d {k, t) - 5M°' d (k, t)] + 

j c d 



E 5 



S ab and X ™>i = h 



(5) 
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where we also used Ajj = Y, a T,b n1A. ab m b . We now as- 
sume that we can choose the moments such that the multi 
relaxation time collision operator is diagonal in moment 
space, i.e. A afc = -5 ab ±. Using T ab (k) = m^e"**"' 1 
and SM° = we thus get the evolution equation of the 
fluctuations in spatial Fourier representation of moment 
space 

5M a (k,t+l) = (12) 

^r ab (fc){(i-l)5M fc (fc,t) + e b (fc,t)}. 

Taking the outer product of SM a with itself, perform- 
ing an ensemble average and substituting r a = 1 - l/r a 
we obtain 

{SM a (k,t + l)SM c (k,t+ 1)) = (13) 
( E E r ab [r b 5M b (k, t) + z b ]r cd 

b d 

[r d 6M d (k,t)+Z d ])- 

For an ideal gas we know the results to be k-independent. 
Henceforth Adhikari et al. only consider the case k = 
at which T ab = S ab . They also invoke stationar- 
ity of equal time correlators {5M a (t + l)SM b (t + 1)) = 
(SM a (t)SM b (t)) and get 

(5M a (t + l)SM c (t + 1)) = r c r a (6M a (t)5M c (t)) + (14) 
r c {SM c (t)C(t)) + r a (8M a {t)i c {t)) + . 

Now, using the fact that the current system state is inde- 
pendent of the noise contribution, i.e. (5M a ^ a ) = 0, they 
obtain 

(£ a £ c ) = (l- r a r c )(SM a SM c ) 

= + » C ~ 1 (6M a 6M c ), (15) 

T a T c 

which acts as the fluctuation dissipation theorem (FDT). 
It relates the noise to the moment fluctuations. What 
is left is finding a prediction for (SM a SM ). Assuming 
the case of the ideal gas [12] they use the fact that the 
distribution functions /,; follow Poisson statistics with a 
mean value and variance of (/j). Thus with Eq. (|5J) they 
get 

(SfiSfj) = (16) 

The back transform to velocity space can now be applied 
to the moment space correlator to obtain 

(5M a 5M b ) = EE< m i(W/i) = (17) 

» j 

This implies that the moment fluctuations and by 
Eq. (|15p the noise terms are generally correlated. How- 
ever, we can decouple these terms by choosing n° = 
m ifi/P because then according to Eq. ([5} 

5»J/?/p=«S o6 (18) 



and thus 

(6M a SM b ) = P 6 ab . (19) 

Of course one has also to show that this is also consis- 
tent with identifying the M a with the hydrodynamic mo- 
ments. For a discussion of this see 

Now that it has been established that the moment fluc- 
tuations can be decoupled according to Eq. (IT9l) we can 
solve Eq. (fTS"]) for the noise amplitude 

e = ±-J p{2 T«-l). (20) 

The actual implementation performes the collision in 
moment space according to Eq. (jB]) where the moments 
M b are constructed at each time step by the standard 
forward transform. The streaming, however, still has to 
happen in velocity space and consequently each update 
involves two matrix transforms. 

Of course, the problem here is that such an orthogo- 
nality condition Eq. (fl~8|) is difficult to fulfill at all times 
and it is not entirely clear which values for p and u we 
have to choose for use in the equilibrium distribution. 
Both Adhikari[2] and Dunweg[3] implicitly assume very 
low flow speeds or the zeroth order expression 

lim/°(p,u)=p Wl , (21) 

thereby avoiding aforementioned problem and simplify- 
ing the orthogonality condition to 

Y l ^m b w t = S ab . (22) 

i 

This implies n" = mfwi and is identical to what is fre- 
quently called the Hermite norm and was originally intro- 
duced by Benzi Q . The orthogonality condition Eq. (|2"2")l 
therefore qualifies the requirements on the transforms in 
addition to the necessity that they preserve hydrodynam- 
ics. An extensive study on the second condition has been 
published in There we found that the Hermite norm 
of Eq. (|22"|) . does indeed also preserve hydrodynamics and 
that, in fact, we are free to add any conserved quantity 
moments to hydrodynamic modes without impacting the 
validity of the hydrodynamic equations. The choice of 
the zeroth order approximation in Eq. (|22p is, however, 
not well documented or motivated in the original litera- 
ture and gives rise to the question whether Galilean in- 
variance violations of the fluctuations result as a conse- 
quence. 



III. GALILEAN INVARIANCE VIOLATIONS IN 
THE HERMITE NORM IMPLEMENTATION 

First we want to evaluate what effect choosing the sim- 
plified norm of Eq. (|2"2"j) has on the Galilean invariance of 
a fluctuating lattice Boltzmann implementation. Here we 



Figure 1: Basis vectors Vi of the D2Q9 scheme used in all 
simulations in this manuscript. 



show the numerical results for an isothermal D2Q9 fluc- 
tuating lattice Boltzmann method with periodic bound- 
ary conditions. Moment space transforms are generated 
with respect to the Hermite norm of Eq. (f2"2"j) . The basis 
vectors Uj are shown in Fig. [TJ All i indices in the fol- 
lowing correspond to these basis vectors. The details of 
the D2Q9 Hermite norm transforms and the equilibrium 
moments are documented in appendix (|A"1) . 

The results in the following were all obtained in a 2D 
lattice Boltzmann simulation of size 21x21. The odd side 
lengths are chosen to avoid the independent conservation 
of momentum components in odd and even lattice sites 
in either dimension. They occur for even side lengths 
because collisions conserve momentum and streaming of 
the densities that constitute momentum and could in- 
teract always moves two lattice sites at once. Conse- 
quently momenta in odd and even numbered lattice sites 
would never interact. We use a large average density of 
po = 10 6 to avoid stability issues due to local negative 
density events. These can occur when the noise £j on 
the distribution functions /, exceeds the value of these 
distribution functions. This is more likely for small p as 
the noise amplitude in moment space Eq. (|20[) is propor- 
tional to sjp. All averages were taken over a simulation 
time of 10 6 iterations after a thermalization phase of 10 5 
iterations to equilibrate the system. 

The fundamental identity that allows us to decouple 
the moment fluctuations is given by Eq. (HHJ). We can 
verify its validity in the simulation directly by measuring 
{SfiSfj) as a function of u x and comparing it to /? and 
Wi of Eq. (flU)) and Eq. (|23p . If the ideal gas hypothesis 
were to hold we would expect Eq. ([To]) to be fulfilled 
independently of u. However, using only the Hermite 
norm Eq. ([22| suggests that we might only find Eq. (fT6|) 
fulfilled to zeroth order, i.e. to the weight factors Wi. 

In Figs. [21 [21 IH we show the simulation results of all 
unique (SfiSfi) correlators as functions of u X: q. We find 
that with increasing velocity u x ^ we do indeed deviate 
strongly from both, the weights Wi, and the equilibrium 
distributions In this implementation the correlators 
approach neither the w; nor the f® and in some cases not 
even an intermediate value. For correlators correspond- 
ing to base velocities without an ^-component ((S/q), 
(<5/|)) the trend opposes that of the f®. In these 



Figure 2: ((<5/o) > in a 21 x 21 D2Q9 fluctuating LB simu- 
lation employing the Hermite norm. We plot Wi and /,°for 
comparison. 
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Figure 3: (((5/0 ) for i = 1...3 in a 21 x 21 D2Q9 fluctuating 
LB simulation employing the Hermite norm. We plot Wi and 
fi for comparison. ((<5/4) 2 ) is not shown as it is identical to 
((<5/2) 2 ) for symmetry reasons. 
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Figure 4: ((Sfij 2 ) for i = 5. ..8 in a 21 x 21 D2Q9 fluctuating 
LB simulation employing the Hermite norm. We plot Wi and 
fi for comparison. ((Sfs) 2 ) and ((Sf?) 2 ) are not shown as 



they appears identical to {(Sfs) 2 ) 
in the scale of this plot. 



and ((Sfs) ) respectively 
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Figure 5: Off-diagonal correlators (SfoSfi) for i = 1...8 in a 
21 x 21 D2Q9 fluctuating LB simulation employing the Her- 
mite norm. {SfoSfi), (SfoSfi), and {SfoSfs) are omitted as 
they behave identical to (<5/o<5/2), (SfoSfo), and {SfoSfs) re- 
spectively. 
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Figure 6: Correlators calculated in the Hermite norm 
{8M a SM a ) normalized to p according to Eq. {I9j in a 21 x 
21 D2Q9 fluctuating LB simulation employing the Hermite 
norm. 
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Figure 7: Linear and quadratic coefficient I and q of all 81 (45 
unique) correlators as a result of fitting {5M a 5M b ){u x ,o) -S ab 
to lu x fi + gu Xi0 . Brighter color indicates larger coefficients. 
Moments were reordered to visually identify correlations bet- 
ter. To accommodate for symbol size the stress moments were 



simplified: IT* 



Tlx 



). The coeffi- 



cient at position (0, 1) in image (a) would correspond to linear 
portion of the {Sj x Sq x ) correlator. Coefficients were measured 
on a 21 x 21 D2Q9 simulation employing the Hermite norm. 
Fit range used was -0.25 <= u x <= 0.25. 



plots and all similar figures in this paper the statisti- 
cal error bars are omitted in the graphs when they are 
smaller than the symbol size. 

In previous publications @, 0| the fluctuations were 
characterized by the fluctuations of the hydrodynamics 
and ghost moments. The corresponding moment corre- 
lators follow directly from the distribution function devi- 
ations according to 

(SM a SM b ) = £ mtm^SfiSfj). (23) 

ij 

and are arguably of more practical importance since they 
represent the fluctuations of the hydrodynamic fields. 

These correlators were expected, in the theory of [||-[(| 
to obey (SM a SM ) = p5 a b- However, for this to work we 
would need (SfiSfj) = Wi in Eq. ([23]) . which is not the 
case for non-zero velocities, as we have shown above. We 
show the observed deviations for the diagonal correlators 
in Fig. [51 Here the correlator of the current in cc-direction, 



(5j x 5j x ), exhibits the largest deviations. 

Note that, while most /j are not symmetric with regard 
to the u x ,q -s- -u Xt o inversion, all the moments are con- 
structed to be either symmetric or antisymmetric under 

U x fi -*■ ~U x fi. 

To obtain some quantitative measure of the depen- 
dency of all 81 (45 unique) correlators in Eq. (|23l) we fit a 
second order polynomial Iu x ,q + qu x to (SM a 6M b )/po - 
S ab . The resulting coefficients I for odd combinations 
and q for even combinations give a rough estimate of the 
deviation of the particular moment correlators and are 
depicted in Fig. [7] We notice in Fig. [7lb) that while 
the quadratic dependency of the correlations on the ve- 
locity is present in several correlators, it is particularly 
apparent on the square correlators. The linear depen- 
dency only appears in cross-correlators which are anti- 
symmetric under u Xt o -*■ -u Xj o as seen in Fig. [7l[a). 

The ensemble averages of the correlation functions 
shown so far do not resolve the length scale dependency 
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of the deviations we observed. To gain some understand- 
ing here we measure the static structure factor 



S^p) = -(Sp(k)S P (--k)). 

Po 



(24) 



the j x momentum correlator 



S k (j x ) = - (5j x (k)S Jx (-k)) , 
Po 



(25) 



at chosen velocities and the momentum cross correlator 



Rk(j*,3v) = - W*(Wv(-k)> (26) 
Po 

at imposed average system velocities u x ,o = 0.0, u Xj q = 
0.1, and u Xt Q = 0.2. We chose Rk(j x ,j y ) in reference 
to Donev et al.'s investigation of the accuracy of finite 
volume schemes (l3j . 

Here Sp(k) = Z x [p( x ) ~ PoK ikx and 5j x (k) = 
E x [ia;( x ) ~ jx,o]e~ lkx are the discrete spatial Fourier 
transforms and £ x i s understood to be the summation 
over all discrete lattice sites. 

In Figs. [H [H and [10] we observe that the correlators 
lose the relatively good agreement with the isotropy re- 
quirement of the ideal gas, i.e. the wave number inde- 
pendence as we increase the velocity. They are sensitive 
to increased velocities and isotropy at the correlations 
is destroyed. Errors are not limited to large k and im- 
pinge on the hydrodynamic (k small) region. Different 
correlators violate isotropy at different length scales and 
directions but we can generalize that the violations for 
certain length scales and spatial directions exceed those 
observed on the level of the ensemble averaged correla- 
tions discussed so far. As an example the density cor- 
relator Su(p) deviates by more than 20% on all length 
scales in the x direction at u x $ = 0.2 in Fig. [5](c) while 
the ensemble average finds a deviation of about 6% in 
Fig. [6] Comparing Figs. H H and [TO] at u^o = 0.2 with 
u x s) = 0.1 we observe that the structure of the anisotropy 
is largely independent of the average system speed al- 
though there are small deviations. Another observation 
is that although (j x j y ) is small compared to other cross 
correlators in Fig. [7] this is mostly due to a fortuitious 
cancellation of errors for different values of k. The abso- 
lute deviations for the (Sj x (k)Sj y (k)) are of similar mag- 
nitude compared to {Sj x (k)Sj x (k)). 

In summary we can clearly see that as function of the 
fluid velocity we observe strong deviations from the iden- 
tities in Eq. (fT^|) and Eq. (JTHJ) and the appearance of off- 
diagonal correlations which are not present in the case of 
u = 0. We conclude that Galilean invariance is indeed 
violated and that the fluctuation-dissipation theorem of 
Eq. (|15p is not longer diagonalized by the simple choice 
of f?/p*Wi inEq. {3). 
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Figure 8: Static structure factor <Sk(p) at different velocities 
measured for the Hermite norm. 



IV. LOCAL VELOCITY DEPENDENT 
TRANSFORMS 



The question now is whether we can alleviate the diffi- 
culties we have encountered by avoiding the approxima- 
tion of /j°(u = 0) = pwi in the normalization condition. 
Removing the velocity dependence in the normalization 
condition could very likely be the source of the Galilean 
invariance violations observed. Instead of using Eq. (|22[) 
we now include the velocity dependence of the equilib- 
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(a) u x ,o = 0.0 




(b) u x ,o = 0.1 




(c) u x ,o = 0.2 



Figure 9: Static structure factor Sk(j x ) at different velocities 
measured for the Hermite norm. 



Figure 10: Cross correlator Rk(jx,j y ) at different velocities 
measured for the Hermite norm. 



rium distribution in Eq. ([18]). The orthogonalization con- orthogonalization with respect to the new scalar product 
dition then becomes 



£<(u)m t 6 (u 



1 1 , x2 1 
1 H U.Vi H (U.Vi) U.U 

6 26> 2 v y 26» 



gab 



(28) 



(27) 

The iterative procedure then follows 

where the velocity u(r,i) is understood to be local to 
the lattice site r. We obtain a new set of transformation 

n—\ 

matrices mf by starting with the physical moments, p, 
jx, j y , Tlxx-yy, Hbj/j ^-xx+yy and perform a Gram-Schmidt 



a ~ & ~ b r0 a 



(29) 



s 




Figure 11: /?(« ,o,u y ) = for all i in the case of the D2Q9 
model. In the area inside the curves fi> for all i. Outside at 
least one /; < and consequently the orthogonalization does 
not find a solution. 



with an intermediate normalization step 



(30) 



With these new matrix elements m° we can define the 
physically relevant moments 



(31) 



One useful side effect of this transform is that the equi- 
librium values for all moments other than the density 
vanish such that 



M a 



p if a = 
otherwise 



(32) 



This is a direct consequence of condition Eq. (f27|) if we 
recognize that M a '° = £im"/°m° = pS a0 because the 
density mode is still the one vector mq = m° = lj. This 
new process does not alter the hydrodynamic limit of 
the lattice Boltzmann method because we only alter the 
moments multiples of u(r) with the conserved quantity 
eigenvectors of the density lj and momentum modes Vi a . 
If we interpret the local velocity u(r) as an arbitrary 
constant we do not alter the hydrodynamic equations at 
all by virtue of our discussion in [9|. We will refer to 
Eq. (|2"Tj) simply as the "/-norm" in the following. 

In order to maintain positive-definiteness of the scalar 
product Eq. (|28p we must be mindful here of the fact 
that the normalization constant needs to be positive at all 
times. The second order expansion of the equilibrium dis- 
tribution Eq. ^ we use here, however, is not. For large 
enough |u| the fifi(p,u,8) < and the orthogonalisation 
has no solution. In Fig.Qj]we show the 0-transition of the 
second order expansion of the equilibrium distribution in 
the case of the D2Q9 model as a function of u. This 
plot shows the accessible velocity range. As long as our 
velocities do not fall outside the central area of Fig. [TT] 



the transformation matrix is guaranteed to be positive 
definite and the Gram-Schmidt will provide a solution. 

The matrix elements m"(u(r)) we obtain are now 
functions of the local velocity u(r) at lattice site r = 
(x,y) . In principle they have to be evaluated at every 
lattice site during every update cycle. We have imple- 
mented a fluctuating LB simulation with these matrices 
and the results are encouraging in that Galilean invari- 
ance violations are significantly smaller. Some results of 
these are shown in Figs. [T^l [T31 and [TU However, even 
in the relatively simple D2Q9 model the matrix elements 
of higher order moments are polynomials of 0(u 16 ) and 
therefore the local evaluation of these matrix elements 
becomes prohibitively costly. Our test implementation 
used between 95% and 99% of the computation time of 
an update cycle in the evaluation of the local transforms. 

One might think that going to the full second order ex- 
pansion of ff might not be necessary and going only to 
first order in u would make the structure of the matrix 
elements significantly simpler. However, working with 
only the first order expansion introduces anisotropy ef- 
fects between the different spatial axis. Removing these 
effectively makes the expressions for the m° even more 
complicated than the regular second order expressions 
where our Gram-Schmidt orthogonalization renders the 
moments isotropic. 

It is, however, not strictly necessary to calculate the 
transforms to machine precision. Judging from our ob- 
servations of the Hermite norm implementation it is suf- 
ficient to calculate tables of the matrix elements on a ve- 
locity grid with velocities u g (g a ) where g a is the grid po- 
sition and use these matrix elements from a look up table 
in the transforms. The benefit is practicality, the pay off 
is that we may not quite obtain the same amount of im- 
provement we might expect to find otherwise. One caveat 
is that we lose the convenient form of the equilibrium mo- 
ments in Eq. (|32[) . In fact the projection of the moments 
in the representation of current local velocity to that of 
the nearest look up table velocity becomes algebraically 
similarly complex as the calculation of the matrix ele- 
ments themselves. However, as we are concerned with 
a second order theory here we choose to only use terms 
of up to O (ug). While we do not change the conserved 
quantities we do change the stress and ghost moments 
at orders O (u 4 ) and higher and thus introduce small er- 
rors. An example of these equilibrium moments and the 
matrix transform elements for D2Q9 can be found in 

The velocity grid spacing for the look up table can be 
relatively coarse. It is helpful if the entire look up ta- 
ble of velocities can fit into the second level cache of the 
CPU the simulation is run on. In our D2Q9 test case 
we typically use a 51 x 51 grid with -0.5 < u g _ x < 0.5, 
-0.5 < Ug^y < 0.5, and Au g = 0.02. Comparing this veloc- 
ity range with Fig. [11] we notice that the corners of this 
square in velocity space falls outside the valid ff (u) > 
range. The matrix elements here are simply evaluated to 
"not a number" and the simulation fails once any one of 
these velocities are reached. In principle one could also 
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Figure 12: {(<5/o) ) in a 21 x 21 D2Q9 fluctuating LB simula- 
tion employing the /-norm with look up tables. Equilibrium 
moments are calculated to third order. {Sfo8fo)f are data- 
points taken from a fully local implementation that foregoes 
the look up table solution. We plot the equilibrium distribu- 
tion /o and the Hermite norm correlator {<5/o5/o}// for com- 
parison. 
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Figure 13: {(Sfi) ) for z = 1...3 in a 21 x 21 D2Q9 fluctuating 
LB simulation employing the /-norm. We plot ff for compar- 
ison. ((<5/4.) 2 } is not shown as it appears identical to {(Sf2) 2 ) 
within the scale of this plot. 



catch outliers in the velocity and just choose the matrix 
elements for a smaller velocity. The moment projection 
would still function. However, this would alter the al- 
gorithm and the results would not be reliable represen- 
tations of the method discussed here. For applications, 
especially at high velocities and low densities it will be 
necessary to include such an exception handling routine. 

One could argue that we might as well have just calcu- 
lated the matrix elements to a lower order directly, forego 
the matrix element look up tables and use the original 
simple equilibrium moments. However, in that case we 
would violate conservation laws and the calculation of the 
2q 2 matrix element polynomials is still significantly more 
expensive than the evaluation of q — d — 1 non-conserved 
moments in a DdQq lattice Boltzmann configuration. 

To evaluate the implementation of the /-norm we per- 
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Figure 14: ((Sfi) ) for i = 5. ..8 in a 21 x 21 D2Q9 fluctuating 
LB simulation employing the /-norm. We plot f° for com- 
parison. {(Sfs) 2 ) and ((<5/t) 2 ) are not shown as they appears 
identical to ((5/s) 2 ) and ((Sfe) 2 } respectively in the scale of 
this plot. 
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Figure 15: (SfoSfi) for i = 1...8 in a 21 x 21 D2Q9 fluctuating 
LB simulation employing the /-norm. 



form the same measurements we did for the Hermite 
norm. We use a D2Q9 ideal gas simulation with peri- 
odic boundaries, and a side length of 21. In Fig. [T2l we 
observe the same (SfoSfo) correlator we did in Fig. [5] 
We find that with the /-norm the trend actually does 
follow the /q prediction and within -0.2 < u x ,o ^ 0.2 we 
are in good agreement with /q but at larger speeds we 
find smaller but noticeable deviations. In Figs. [TSJ Q~l] 
we find much better agreement for all other distribution 
function correlation functions for the /-norm compared 
to the Hermite norm in Figs. [31 andUJ Again we notice 
very good agreement for \u x \ < 0.2. 

The remaining deviations from the equilibrium distri- 
butions we find with the /-norm are not an artifact of 
either the look up table method or the third order ex- 
pansion of the equilibrium moments. We performed the 
same measurement with the fully locally orthogonalized 
set of transforms, albeit with fewer data points due to 
the much higher computational effort involved. (SfiSfi) t 
in Figs. [TSJ [TSJ and [141 indicate that the deviations from 
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Figure 16: Correlators {5M a 8M a ) / normalized to p accord- 
ing to Eq. (fT§)) in a 21 x 21 D2Q9 fluctuating LB simulation 
employing the /-norm. 
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Figure 17: Correlators {SM a 8M a ) normalized to p accord- 
ing to Eq. (fT9)l measured in a 21 x 21 D2Q9 fluctuating LB 
simulation employing the Hermite norm. 



the equilibrium distributions can indeed not be explained 
with either the look up table method or the cut off on 
the equilibrium moments as the results obtained form 
the look up table method with third order equilibrium 
moments appears to be consistent from the fully locally 
orthogonalized /-norm. 

Measuring the moment space correlators in the /-norm 
poses an interesting question. Do we measure with re- 
spect to the Hermite norm or the /-norm and in the case 
of the latter with respect to which velocity? To answer 
this question we conduct a thought experiment. 8M a 
should be Galilean invariant for any a, in particular the 
momentum components. In the Hermite norm we have 

Sjx = Y, m ifi ~ zL m ifi = ^ (pu x - p Q U xfi ) (33) 



we have Sj x = Sj x = \/3pu x . Introducing a constant ve- 
locity offset -uo should leave Sj x Galilean invariant, i.e. 
we expect u -* u - uo ■ If we now interpret Uo as such 
an offset the Hermite norm is clearly not Galilean in- 
variant under velocity offsets as it introduces an extra 
u x ,o (po ~ p) m Eq. (|33|) whereas the /-norm in Eq. (|34[) 
behaves as required. Consequently we use the /-norm as 
it provides the correct measurements that leave the SM a 
invariant under Galilean transformations. Furthermore 
we measure with respect to the average system velocity 
Uo and average density po. Measuring with respect to 
the local velocity u and density p is nonsensical as 6p = 
and (5j = in this case. We thus use the /-norm such 
that fh^fh\{fi) = 5 ab where we make the approximation 
of Eq. ©(/<> = /P(po,uo). 

Much like the distribution function correlators the 
moment correlators ((SAl a ) 2 } shown in Fig. [TBI exhibit 
significant improvement compared to those of the Her- 
mite norm in Fig. [BJ This improvement is smaller than 
the general trend of the distribution function correla- 
tors would imply for some modes. In particular the 
((Sp) 2 ), ((6Yl xx - yy ) 2 ), and ((Sj y ) 2 ) correlators deviate 
significantly for larger u x . Their overall decrease is about 
1/3 compared to the Hermite norm. To make a valid com- 
parison between moment correlators computed in the /- 
norm and the Hermite norm one needs to ensure that 
for both measurements the moments are obtained in the 
same way. We therefore measure the moments obtained 
in a Hermite norm simulation with the /-norm evaluated 
at Uo in Fig. [T71 We observe that for all moments but 
(SpSp) and {Sj y Sj y } the deviations are larger than those 
measured in the Hermite norm. 

Linear and quadratic fit coefficients for all moment cor- 
relators (SM a SM b ) in Fig. ITS1 show significant improve- 
ment as well. We notice that in particular the coefficients 
I that apply to those off-diagonal correlators that have a 
linear dependence on u x at least a factor of 13 smaller 
than those measured in the Hermite norm case shown in 
Fig. [7] (a). We also observe a decrease of the quadratic 
term q but in line with the observations of Fig.[Tj)]the co- 
efficients corresponding to some correlators decrease less 
compared to the ones observed in the Hermite norm in 
Fig. E| (b): {(5p) 2 ) from 1.9 to 0.47, {(SU xx ^ yy ) 2 ) from 
1.6 to 0.54, and {(Sj y ) 2 ) from 1.14 to 0.75. 

These findings are confirmed by the structure factor 
plots for the /-norm in Figs. [1]JJ [201 and which for 
non- vanishing fixed velocity u x ,o & re significantly smaller 
than the one measured for the Hermite norm at the same 
velocity in Figs. GO El and OS 



and for the /-norm 

5% = E mffi - £ mifi = V3p (u x - u xfi ) . (34) 

i i 

Again u is the mean velocity in the system and u the 
local velocity at a given lattice site. If we set Uq = 



We can conclude that employing the /-norm signifi- 
cantly reduces the Galilean invariance effects observed on 
the Hermite norm implementation. The look up tables 
provide a practically feasible approach to implementing 
the /-norm at a performance loss of about 20 %. All the 
measurements here were performed on a single CPU. 
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Figure 18: Linear and quadratic coefficient I and q of all 81 (45 
unique) correlators as a result of fitting (SM a 5M b )(u x ,o)-5 ab 
to lu x fi + qu x> Q. Brighter color indicates larger coefficients. 
Moments were reordered to visually identify correlations bet- 
ter. To accommodate for symbol size the stress moments were 



simplified: LT X 



,n + =a x 



The coeffi- 



cient at position (0, 1) in image (a) would correspond to linear 
portion of the (Sj x 5q x ) correlator. Coefficients were measured 
on a 21 x 21 D2Q9 simulation employing the /-norm with look 
up tables, u g = 0.02. Fit range used was -0.25 <= u x <= 0.25. 



V. CONCLUSION AND OUTLOOK 

The current standard implementation of thermal fluc- 
tuations in an isothermal ideal gas was tested for Galilean 
invariance violations. We found that with non zero av- 
erage velocity the moment space covariance matrix of 
Eq. (fTT)|) is neither diagonal nor are the diagonal elements 
unity as predicted and required by the derivation of the 
FDT in both @| and 0]. We identified an approxima- 
tion in the orthogonality condition that defines the mo- 
ment space transforms Eq. (|18[) as the likely source of 
the Galilean invariance violations as it directly removes 
an otherwise necessary velocity dependence from the mo- 
ment space transforms. The approximation allows for the 
use of Hermite norm to define the moment space trans- 
forms. However, to recover Galilean invariance at least 
to some degree requires the matrix transforms to be lo- 
cally velocity dependant, i.e. unique to every lattice site 
and the Hermite norm is no longer applicable. This led 




Figure 19: Static structure factor Sk(p) at different velocities 
measured for the /-norm with the look up table and Au g = 
0.02. 



us to introduce a novel variant of the lattice Boltzmann 
method. We find that using the local fully velocity depen- 
dent /-norm to machine precision in a straight forward 
manner to be computationally impractical. Evaluating 
the individual matrix elements leads to an overhead in 
computational cost of > 2000% in evaluating the individ- 
ual matrix elements. However, as the Galilean invariance 
violations scale quadratically for most moments it is fea- 
sible to generate look up tables for the matrix elements 
on a velocity grid. This requires to projection of the 
equilibrium moments into the look up table reference ve- 
locity. This look up table approach provides comparable 
benefits to the locally orthogonalized transforms but at 
only a 20% loss of computation time. All the simula- 
tions presented here were performed in a example D2Q9 
implementation. However, all calculations and considera- 
tions discussed can easily be generalized to other models. 
We provide a Mathematica notebook [14( that contains 
the necessary calculations done for the D2Q9 model used 
here. This new method is poentially important for non- 
equilibrium situations when locally varying flow fields 
exist which is the standard realm of lattice Boltzmann 
simulations. 
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Appendix A: Hermite norm D2Q9 



For D2Q9 the equilibrium distribution employed is 
*iven by Eq. © with 9=1/3 



/°(P,M) = P w i 



9 2 3 
1 + 3u.Vi + — (u.Vi) - -u.u 



(Al) 



The weights are given by 



I if * = 

| if z = 1,2,3,4 

t^t if i = 5, 6, 7, 8 



(A2) 



Figure 20: Static structure factor Sk(jx) at different velocities 
measured for the /-norm with the look up table and Au g = 
0.02. 



In the case of the simple Hermite norm Eq. (|22|) it is fea- 
sible to show the transformation matrices. The forward 
transform reads 
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Figure 21: Cross correlator Rk(jx,j y ) at different velocities 
measured for the /-norm with the look up table and Au g = 
0.02. 



Likewise the back transform from moment space to velocity space is given by 
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where n" = mfwi The corresponding equilibrium moments M a ' are obtained directly by applying the forward trans- 
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form to the equilibrium distribution. In the Hermite norm we find 
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